Code
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mntFigure 4 and supplementary figure 4
1 Intro
We performed histone PTMs DIA MS and I analysed the raw files with EpiProfile v2.1 (Yuan et al. 2018). Here I parse the output and analyse the PTMs values with DEP.
Acknowledgement: I’d like Zuo-Fei Yuan and Simone Sidoli from Ben Garcia lab for their help with data analysis and sample preparation respectively.
2 Set Up
To fetch files on the CRG cluster I mount the server on my local computer with the following command in the terminal.
2.1 Packages
Load packages required for the analysis and suppress any message. Check the Section 5 section at the end for package versions.
Code
library(dplyr, warn.conflicts = F, quietly = T)
library(readr)
library(readxl)
library(DESeq2)
library(csaw)
library(pheatmap)
library(DEP)
library(Cairo)
# library(knitr)
library(DT)Load my R package niar to use my custom made function to parse vast-tools output.
Code
if ( ('niar' %in% .packages(all.available = TRUE)) == TRUE ) {
library('niar')
} else if ( ('niar' %in% .packages(all.available = TRUE)) == FALSE ){
message("The package niar is not installed so I'm trying to load it from ",
"the local repository of the package")
if ( dir.exists('~/mnt/narecco/software/R/niar' ) ) {
devtools::load_all(path = '~/mnt/narecco/software/R/niar')
} else {
stop("Can't find the local repo of the niar package! ",
"You must install it with:\n",
"devtools::install_github('Ni-Ar/niar') ")
}
} else{
stop("Can't understand if 'niar' package was installed beforehand")
}The package niar is not installed so I'm trying to load it from the local repository of the packageℹ Loading niar2.2 Functions
Define a plot style.
Code
theme_classic(base_size = 6.3, base_family = "Arial") +
theme(plot.title = element_text(size = 8, vjust = -1, hjust = 0),
plot.background = element_blank(),
strip.background = element_blank(),
panel.background = element_blank(),
axis.line = element_line(linewidth = 0.1),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 5),
axis.text = element_text(colour = "black"),
axis.text.x = element_text(margin = margin(t = 0, unit = "mm")),
axis.ticks.x = element_blank(),
axis.ticks.y = element_line(linewidth = 0.1),
axis.ticks.length.y = unit(1, "mm"),
legend.margin = margin(l = 0, t = 0, unit = "mm"),
legend.box.margin = margin(l = 0, unit = "mm"),
legend.box.background = element_blank(),
legend.box.spacing = unit(0.5, "mm"),
# legend.title = element_blank(),
legend.text = element_text(margin = margin(l = -10, unit = "mm")),
legend.text.align = 0,
legend.key.height = unit(2, "mm"),
legend.key.width = unit(7.5, "mm"),
legend.key = element_blank(),
legend.background = element_blank(),
legend.spacing.x = unit(3, 'mm')
) -> alluvial_thm I modify the geom_stratum() function just to add a very small linewidth the the black boders of the colums. With linewidth = .05, I don’t have to change them manually later
Small function to parse the PTM abundance as 10 to the power scientific notion.
Code
scientific_10 <- function(x) {
parse(text=gsub("e\\+", " %*%10^", scales::scientific_format()(x)))
}2.3 Directories & File Paths
Here I organise all the file and directories paths I need to run the analysis and define where to save the processed tables and figures.
Code
oneDrive_Dir <- file.path("~/OneDrive - CRG - Centre de Regulacio Genomica/Suz12_AS_project")
code_dir_fig <- file.path(oneDrive_Dir, "_Code/Fig4")
tbl_dir_fig4 <- file.path(code_dir_fig, "tables")
pdf_dir_fig4 <- file.path(code_dir_fig, "pdfs")
if (!dir.exists(pdf_dir_fig4)) { dir.create(pdf_dir_fig4, recursive = T) }
if (!dir.exists(tbl_dir_fig4)) { dir.create(tbl_dir_fig4, recursive = T) }Code
metadata_path <- file.path(tbl_dir_fig4, 'Histone_PTMs_metadata.xlsx')
WT_Dex4_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet1_Run2_DIA_histone_ratios_Suz12dExon4.xls",
full.names = T)
Resce_KO_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet2_Run1_DIA_histone_ratios_Suz12KO.xls",
full.names = T)
stopifnot(file.exists(metadata_path))
stopifnot(file.exists(WT_Dex4_histone_ratio_path))
stopifnot(file.exists(Resce_KO_histone_ratio_path))3 Main Figure Panels
Import metadata.
Code
mtdt <- read_excel(metadata_path) |>
setNames(c("Sample", "Short_Name", "Clone_Name", "Condition", "Replicate"))
mtdt$Condition <- factor(mtdt$Condition,
levels = c("WT", "CSexon4", "Dexon4", "KO", "KO_L", "KO_S") )3.1 H3 K27-K40 peptide modification ratios
Import DIA data quantification using automatic EpiProfile v2.1 quantifications.
Code
read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set1
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set2
data_dia_auto <- rbind(dia_auto_set1, dia_auto_set2)Tidy up the information for histone peptide H3 K27- K36. Here I split the H3 K27 and K36 signal coming form the same peptide. In this way I can pile up the all methyl-forms.
Code
data_dia_auto |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(K27 = case_when( is.na(Sub_PTM1) & grepl(pattern = "^K27", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
mutate(K36 = case_when( is.na(Sub_PTM2) & grepl(pattern = "^K36", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
relocate(K27, K36, .after = Sub_PTM2) |>
unite(col = "K27", c("Sub_PTM1", "K27"), sep = "", na.rm = T, remove = T ) |>
unite(col = "K36", c("Sub_PTM2", "K36"), sep = "", na.rm = T, remove = T ) |>
mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
# When one mark is not found I label it unmodified ('unmod') but it could also
# be labelled as "only the other mark".
# mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "Only\nK36me",
# grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
# mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "Only\nK27ac-me",
# grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
relocate(K27, K36, .after = PTM) |>
mutate(Histone = factor(Histone, levels = c("H3", "H3.3") ) ) -> data_K27_36I can also deconvolve the K27 and K36 signals from the same peptide and pile-up each signal.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) )|>
select(Histone, Condition, Average_Ratio, K27, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K27, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K27) |> ungroup() -> data_K27Same as above but for K36.
Code
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
select(Histone, PTM, Condition, Average_Ratio, K36, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K36, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K36) |> ungroup() -> data_K36For specific histone PTMs one can explore the results in the following tables.
Overview of H3 and H3.3 K27 acetylation and methylation quantifications
Code
datatable(data_K27, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )Overview of H3 and H3.3 K36 methylation quantifications
Code
datatable(data_K36, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )3.1.1 H3 K27
Plot figure 4B panel with the pileup of figure H3K27 methylation and acetylation.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |> # 3 replicates
subset(Histone %in% "H3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_K27
p_Alluvial_K27 
Save to pdf.
Code
ggsave(filename = "Fig4B_H3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)3.2 Differential abbundance analysis of histone PTMs with DEP
Use the area under the peak as counts to be analysed using DEP framework.
Code
rbind( read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path),
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) ) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) -> data_dia_auto_not_tidyOne small caveat is that the areas under the peaks are very large numbers or at least larger than 2147483647 which is the highest number R can store as integer. So to convert the areas to integers from numeric type one could use the package bit64 that stores integers in 64 bits instead of the default 32 bit. My solution was however simpler, I divided the all table by 1000 to scale down by 3 order of magnitude. I then convert the dataframe to matrix.
I tried to add a pseudocounts to the data but it actually doesn’t help with the analysis and is better to input the areas as they are. The function make_unique is the first step of the analysis of DEP and creates a unique data frame.
Code
data_dia_auto_not_tidy |>
tidy_hPTMs() |>
select(First_Col, PTM, Sample, Area) |>
# reduce all areas by a 1000 to avoid having to coerc to integer big numbers
mutate(Area = Area / 1e4) |>
pivot_wider(id_cols = c(First_Col, PTM), names_from = Sample, values_from = Area) |>
# Coerce the areas to integers
mutate(across(.cols = ends_with(c("_1", "_2", "_3")), .fns = as.integer )) |>
make_unique(names = "First_Col", ids = "PTM") -> data_uniqueCreate a metadata compatible for DEP set the KO as the reference sample as I’m interesting in the rescue efficiency relative to the KO.
Code
# Set metadata colnames to be as `make_se` wants them to be
colnames(mtdt) <- c("Sample", "label", "Clone_Name", "condition", "replicate")
mutate(mtdt, condition = factor(condition, levels = c("WT", "CSexon4", "Dexon4",
"KO", "KO_L", "KO_S"))) |>
print(n = 1 ) -> mtdt_dep# A tibble: 18 × 5
Sample label Clone_Name condition replicate
<chr> <chr> <chr> <fct> <dbl>
1 2022MZ006_IVMO_001_02_25pto_DIA WT_1 ZWT WT 1
# ℹ 17 more rowsFind the columns that contain the samples areas in the data_unique dataframe and create a SummarizedExperiment object to be used in for the DEP analysis. The first step is to do a log2 transformation of the area counts.
Code
Samples_area_columns <- grep(pattern = paste(mtdt_dep$label, collapse = "|"), colnames(data_unique))
data_se <- make_se(data_unique, columns = Samples_area_columns, expdesign = mtdt_dep)Now the data_se contains the modifications in the rows (assuming they are proteins) and in the column the area.
Code
data_filt <- filter_missval(data_se, thr = 0)Check PTMs frequency across samples.
Code
plot_frequency(data_filt)
Overall most modifications are found in all samples.
Code
data_norm <- normalize_vsn(data_filt)Verify the fit.
Code
meanSdPlot(x = data_norm)
Seems okay.
Check normalisation effect on data.
Code
plot_normalization(data_se, data_filt, data_norm) +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1")) +
theme(axis.text = element_text(color = 'black'),
plot.background = element_blank(),
panel.background = element_blank())
Check missing values.
Code
plot_missval(data_norm)
I don’t see obvious missing values trends, besides only that sample 3 of CSexon 3 has less PTMs detected. I don’t think it’ necessary to impute missing values.
Code
# Test every sample versus control
data_diff <- test_diff(data_norm, type = "control", control = "WT")Warning in test_diff(data_norm, type = "control", control = "WT"): Missing
values in 'data_norm'Tested contrasts: CSexon4_vs_WT, Dexon4_vs_WT, KO_vs_WT, KO_L_vs_WT, KO_S_vs_WTWarning: Partial NA coefficients for 5 probe(s)Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!
Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!
Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!
Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!
Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!Significance thresholds:
P-value <= 0.05
|Log2 Fold Change| >= 1.5
Code
dep <- add_rejections(data_diff, alpha = 0.05, lfc = 1.5 )Check sample clustering a Gower’s distance matrix.s
Code
plot_dist(dep, significant = TRUE, pal = "PuOr") 
Check sample clustering.
Code
plot_heatmap(dep, type = "centered", kmeans = F, clustering_distance = "gower",
col_limit = 4, show_row_names = TRUE,
indicate = c("condition", "replicate"))Warning: Missing values in 'dep'. Using clustering_distance = 'gower'
Get results.
Code
data_results <- get_results(dep)
data_results |>
as_tibble() |>
select(-significant) |>
pivot_longer(cols = !c("name", "ID"),
names_to = c("contrast", "Feature"),
names_sep = c("_(vs_WT_|centered)"),
values_to = "Value"
) |>
mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
pivot_wider(names_from = "Feature", values_from = "Value") |>
dplyr::select(-ID) -> res_all_tidy
colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'Extract WT sample area of each PTM, to be used for the volcano plot.
Code
data_dia_auto |>
subset(Condition == "WT") |>
group_by(First_Col) |>
mutate(WT_Mean_Area = mean(Area) + 1 ) |>
relocate(WT_Mean_Area, .after = PTM) |>
select(c(First_Col, WT_Mean_Area)) |>
unique() -> WT_areaGet the PTMs info columns.
Code
histone_ptms_anno <- unique(data_dia_auto[, 1:10])Add all the info the the tidy results dataframe.
Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )
res_all |>
mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
arrange(desc(Histone), Peptide_Start, Peptide_End) |>
ungroup() -> res_allExport all tidy results for supplementary table 3 as tab delimited, or as csv with UTF-8 Byte order mark which indicates to Excel the csv is UTF-8 encoded.
Code
write_delim(x = res_all, file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.tab'),
delim = '\t', col_names = T, append = F, quote = 'none', na = "NA",
progress = F, escape = 'none')
write_excel_csv(x = res_all,
file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.csv'),
delim = ',', col_names = T, append = F, quote = 'none',
na = "NA", progress = F, escape = 'none')Interactively explore the table
Code
res_all |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
datatable(rownames = FALSE, filter = 'top',
options = list(pageLength = 3, autoWidth = TRUE) )3.2.1 PTMs centred abundances
Plot centred abundance
Code
res_all |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
subset(!is.na(centered)) |>
ggplot(aes(x = centered, y = PTM, fill = contrast)) +
facet_grid(rows = vars(Modification), scales = 'free' ) +
geom_vline(xintercept = 0, linewidth = 0.4) +
geom_col(width = 0.75,
position = position_dodge(width = 0.75, preserve = 'single')) +
geom_point(aes(size = WT_Mean_Area),
shape = 21, position = position_dodge(width = 0.75),
stroke = 0.2) +
labs(x = "Centered PTM abundance") +
scale_x_continuous(limits = c(-5, 5), oob = scales::oob_squish ) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1"),
name = "Sample") +
theme_classic() +
theme(strip.text = element_blank(),
strip.background = element_blank(),
axis.text = element_text(colour = "black"),
axis.title.y = element_blank(),
legend.key.size = unit(1, 'mm'),
panel.grid.major.y = element_line()) -> p_centered_K27_K36
p_centered_K27_K36
3.2.2 Volcano plot ∆ex4 vs WT
Prepare data for plot.
Code
res_Dex4 <- subset(res_all, contrast == "Dexon4" & !is.na(p.val))
res_Dex4 <- res_Dex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot figure 4C.
Code
ggplot(res_Dex4 ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 5.5, y = 0.10, label = "∆ex4",
colour = "black", fill = "goldenrod", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -5.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.09 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.1 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-6.4, 6.4), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.91, 0.87),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_Dex4
p_Dex4
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "Fig4C_Dex4_vs_WT_Volcano_w_labels.pdf",
plot = p_Dex4, device = cairo_pdf, units = "cm",
height = 5.4, width = 7.0)3.2.3 Other volcano plots
Prepare data of CSex4 vs WT
Code
res_CSex4 <- subset(res_all, contrast == "CSexon4" & !is.na(p.val))
res_CSex4 <- res_CSex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Note: nothing is differentially expressed in this contrast
Code
table(res_CSex4$Direction)
None
159 Prepare data of KO vs WT
Code
res_KO <- subset(res_all, contrast == "KO" & !is.na(p.val))
res_KO <- res_KO |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot KO vs WT
Code
ggplot(res_KO ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 8.5, y = 0.10, label = "KO",
colour = "white", fill = "gray16", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -12.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-15, 10), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = 'right',
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_KO
p_KO
Note: these WT and KO samples were processed at different moment, so there could be some batch effect difference that limits the statistical power.
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "KO_vs_WT_Volcano_w_labels.pdf",
plot = p_KO, device = cairo_pdf, units = "cm",
height = 12, width = 10)
Things I could still add to this analysis
I could still add:
the rescues volcano plots vs the KO and KO+S vs KO+L
rescue efficiency
heatmap with all PTMs.
4 Supplementary Figure Panels
4.1 K27-K40 peptide modification ratios
4.1.1 H3.3K27
There’s no H3.3 signal in the main figure, here I plot only H3.3 signal.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone %in% "H3.3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3.3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_H33K27
p_Alluvial_H33K27 
Save to pdf.
Code
ggsave(filename = "FigS4B_H3.3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_H33K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)4.1.2 H3 and H3.3 K36
Plot the pileup of H3 & H3.3 K36 methylation.
Code
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone == c("H3", "H3.3") ) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K36, alluvium = K36,
y = Ratio_Sum, fill = K36) +
facet_wrap(~ Histone, scales = 'free_y') +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, linewidth = 0.1, colour = "black") +
geom_fit_text(aes(label = K36 ), stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"), size = 7, family = "Arial",
padding.y = grid::unit(0.2, "mm"),
show.legend = F) +
scale_fill_manual(values = c("K36me1" = "#87B7AA", "K36me2" = "#8FA9C2",
"K36me3" = "#8386B7", "unmod" = "gray" ) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO', 'KO+L', 'KO+S')
) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.title = element_blank()) -> p_Alluvial_K36
p_Alluvial_K36
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "FigS4C_H3_H3.3_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K36, device = cairo_pdf, units = "cm",
width = 12.0, height = 5)4.2 Extra plots
Here I plot extra histone PTMs that were not included in the publication.
4.2.1 H3 K4
Plot H3K4 modifications.
Code
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Peptide_Start >= 3 & Peptide_End <= 8) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Demuth", direction = 1, n = 5), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -9, unit = "mm")),
legend.key.width = unit(6.5, "mm") ) -> p_Alluvial_K4
p_Alluvial_K4
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "H3K4_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K4, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)4.2.2 H3 K9
Plot all combinatorial K9 and K14 modifications on histone H3.
Code
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Histone == "H3") |>
subset(Peptide_Start >= 9 & Peptide_End <= 17) |>
mutate(Average_Ratio = mean(Ratio) ) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hokusai1", direction = -1, n = 15), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -12.0, unit = "mm")),
legend.key.width = unit(11.2, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K9
p_Alluvial_K9
Code
ggsave(path = pdf_dir_fig4, filename = "H3K9_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K9, device = cairo_pdf, units = "cm",
width = 6.4, height = 4.85)4.2.3 All combinations of H3/H3.3 K27 K36
Plot all combinatorial K27 and K36 modifications on histone H3 and H3.3 variant.
Code
data_dia_auto |>
group_by(Condition, Modification, Histone) |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
facet_wrap(~ Histone) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hiroshige", direction = 1, n = 15) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S'), name = 'H3 & H3.3') +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -14.5, unit = "mm")),
legend.key.width = unit(13.75, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K27_K36
p_Alluvial_K27_K36
This is the best representative picture of K27-K40 modifications relative abundance.
Code
ggsave(path = pdf_dir_fig4, filename = "H3_H3.3_K27_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K27_K36, device = cairo_pdf, units = "cm",
width = 14, height = 5)End analysis.
5 Session Info
Code
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.1.2 (2021-11-01)
os macOS Catalina 10.15.7
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Madrid
date 2023-09-04
pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
! package * version date (UTC) lib source
ade4 1.7-22 2023-02-06 [1] CRAN (R 4.1.2)
affy 1.72.0 2021-10-26 [1] Bioconductor
affyio 1.64.0 2021-10-26 [1] Bioconductor
annotate 1.72.0 2021-10-26 [1] Bioconductor
AnnotationDbi 1.56.2 2021-11-09 [1] Bioconductor
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
Biobase * 2.54.0 2021-10-26 [1] Bioconductor
BiocFileCache 2.2.1 2022-01-23 [1] Bioconductor
BiocGenerics * 0.40.0 2021-10-26 [1] Bioconductor
BiocManager 1.30.22 2023-08-08 [1] CRAN (R 4.1.2)
BiocParallel 1.28.3 2021-12-09 [1] Bioconductor
biomaRt 2.50.3 2022-02-06 [1] Bioconductor
Biostrings 2.62.0 2021-10-26 [1] Bioconductor
bit 4.0.5 2022-11-15 [1] CRAN (R 4.1.2)
bit64 4.0.5 2020-08-30 [1] CRAN (R 4.1.0)
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.0)
blob 1.2.4 2023-03-17 [1] CRAN (R 4.1.2)
bslib 0.5.1 2023-08-11 [1] CRAN (R 4.1.2)
cachem 1.0.8 2023-05-01 [1] CRAN (R 4.1.2)
Cairo * 1.6-1 2023-08-18 [1] CRAN (R 4.1.2)
callr 3.7.3 2022-11-02 [1] CRAN (R 4.1.2)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
circlize 0.4.15 2022-05-10 [1] CRAN (R 4.1.2)
cli 3.6.1 2023-03-23 [1] CRAN (R 4.1.2)
clue 0.3-64 2023-01-31 [1] CRAN (R 4.1.2)
cluster 2.1.4 2022-08-22 [1] CRAN (R 4.1.2)
codetools 0.2-19 2023-02-01 [1] CRAN (R 4.1.2)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.1.2)
ComplexHeatmap 2.10.0 2021-10-26 [1] Bioconductor
crayon 1.5.2 2022-09-29 [1] CRAN (R 4.1.2)
crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.1.0)
csaw * 1.28.0 2021-10-26 [1] Bioconductor
curl 5.0.2 2023-08-14 [1] CRAN (R 4.1.2)
DBI 1.1.3 2022-06-18 [1] CRAN (R 4.1.2)
dbplyr 2.3.3 2023-07-07 [1] CRAN (R 4.1.2)
DelayedArray 0.20.0 2021-10-26 [1] Bioconductor
DEP * 1.16.0 2021-10-26 [1] Bioconductor
desc 1.4.2 2022-09-08 [1] CRAN (R 4.1.2)
DESeq2 * 1.34.0 2021-10-26 [1] Bioconductor
devtools 2.4.5 2022-10-11 [1] CRAN (R 4.1.2)
digest 0.6.33 2023-07-07 [1] CRAN (R 4.1.2)
doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.1.2)
dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.1.2)
DT * 0.29 2023-08-29 [1] CRAN (R 4.1.2)
edgeR 3.36.0 2021-10-26 [1] Bioconductor
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
evaluate 0.21 2023-05-05 [1] CRAN (R 4.1.2)
fansi 1.0.4 2023-01-22 [1] CRAN (R 4.1.2)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.1.2)
fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.1.2)
fdrtool 1.2.17 2021-11-13 [1] CRAN (R 4.1.0)
filelock 1.0.2 2018-10-05 [1] CRAN (R 4.1.0)
forcats 1.0.0 2023-01-29 [1] CRAN (R 4.1.2)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.1.2)
foreign 0.8-84 2022-12-06 [1] CRAN (R 4.1.2)
fs 1.6.3 2023-07-20 [1] CRAN (R 4.1.2)
genefilter 1.76.0 2021-10-26 [1] Bioconductor
geneplotter 1.72.0 2021-10-26 [1] Bioconductor
generics 0.1.3 2022-07-05 [1] CRAN (R 4.1.2)
GenomeInfoDb * 1.30.1 2022-01-30 [1] Bioconductor
GenomeInfoDbData 1.2.7 2023-08-29 [1] Bioconductor
GenomicRanges * 1.46.1 2021-11-18 [1] Bioconductor
GetoptLong 1.0.5 2020-12-15 [1] CRAN (R 4.1.0)
ggalluvial 0.12.5 2023-02-22 [1] CRAN (R 4.1.2)
ggfittext 0.10.0 2023-04-04 [1] CRAN (R 4.1.2)
ggplot2 3.4.3 2023-08-14 [1] CRAN (R 4.1.2)
ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.1.2)
ggseqlogo 0.1 2017-07-25 [1] CRAN (R 4.1.0)
GlobalOptions 0.1.2 2020-06-10 [1] CRAN (R 4.1.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.2)
gmm 1.8 2023-06-06 [1] CRAN (R 4.1.2)
gtable 0.3.4 2023-08-21 [1] CRAN (R 4.1.2)
hexbin 1.28.3 2023-03-21 [1] CRAN (R 4.1.2)
hms 1.1.3 2023-03-21 [1] CRAN (R 4.1.2)
htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.1.2)
htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.1.2)
httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.1.2)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.1.2)
impute 1.68.0 2021-10-26 [1] Bioconductor
imputeLCMD 2.1 2022-06-10 [1] CRAN (R 4.1.2)
IRanges * 2.28.0 2021-10-26 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.1.2)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.1.2)
KEGGREST 1.34.0 2021-10-26 [1] Bioconductor
knitr 1.43 2023-05-25 [1] CRAN (R 4.1.2)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
later 1.3.1 2023-05-02 [1] CRAN (R 4.1.2)
lattice 0.21-8 2023-04-05 [1] CRAN (R 4.1.2)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.1.2)
limma 3.50.3 2022-04-07 [1] Bioconductor
locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.1.2)
magick 2.7.5 2023-08-07 [1] CRAN (R 4.1.2)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.2)
MALDIquant 1.22.1 2023-03-20 [1] CRAN (R 4.1.2)
MASS 7.3-60 2023-05-04 [1] CRAN (R 4.1.2)
Matrix 1.6-1 2023-08-14 [1] CRAN (R 4.1.2)
MatrixGenerics * 1.6.0 2021-10-26 [1] Bioconductor
matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.1.2)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.1.0)
metapod 1.2.0 2021-10-26 [1] Bioconductor
MetBrewer 0.2.0 2022-03-21 [1] CRAN (R 4.1.2)
mime 0.12 2021-09-28 [1] CRAN (R 4.1.0)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.0)
msa 1.26.0 2021-10-26 [1] Bioconductor
MsCoreUtils 1.6.2 2022-02-24 [1] Bioconductor
MSnbase 2.20.4 2022-01-16 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
mvtnorm 1.2-3 2023-08-25 [1] CRAN (R 4.1.2)
mzID 1.32.0 2021-10-26 [1] Bioconductor
mzR 2.28.0 2021-10-27 [1] Bioconductor
ncdf4 1.21 2023-01-07 [1] CRAN (R 4.1.2)
R niar * 0.3.0 <NA> [?] <NA>
norm 1.0-11.1 2023-06-18 [1] CRAN (R 4.1.2)
patchwork 1.1.3 2023-08-14 [1] CRAN (R 4.1.2)
pcaMethods 1.86.0 2021-10-26 [1] Bioconductor
pheatmap * 1.0.12 2019-01-04 [1] CRAN (R 4.1.0)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.1.2)
pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.1.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
pkgload 1.3.2.1 2023-07-08 [1] CRAN (R 4.1.2)
plyr 1.8.8 2022-11-11 [1] CRAN (R 4.1.2)
png 0.1-8 2022-11-29 [1] CRAN (R 4.1.2)
preprocessCore 1.56.0 2021-10-26 [1] Bioconductor
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.1.0)
processx 3.8.2 2023-06-30 [1] CRAN (R 4.1.2)
profvis 0.3.8 2023-05-02 [1] CRAN (R 4.1.2)
progress 1.2.2 2019-05-16 [1] CRAN (R 4.1.0)
promises 1.2.1 2023-08-10 [1] CRAN (R 4.1.2)
ProtGenerics 1.26.0 2021-10-26 [1] Bioconductor
ps 1.7.5 2023-04-18 [1] CRAN (R 4.1.2)
psychTools 2.3.6 2023-06-18 [1] CRAN (R 4.1.2)
purrr 1.0.2 2023-08-10 [1] CRAN (R 4.1.2)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.1.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.1.2)
Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.1.2)
RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.1.2)
readr * 2.1.4 2023-02-10 [1] CRAN (R 4.1.2)
readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.1.2)
remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.1.2)
rjson 0.2.21 2022-01-09 [1] CRAN (R 4.1.2)
rlang 1.1.1 2023-04-28 [1] CRAN (R 4.1.2)
rmarkdown 2.24 2023-08-14 [1] CRAN (R 4.1.2)
rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.2)
Rsamtools 2.10.0 2021-10-26 [1] Bioconductor
RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.1.2)
rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.1.2)
S4Vectors * 0.32.4 2022-03-29 [1] Bioconductor
sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.1.2)
sass 0.4.7 2023-07-15 [1] CRAN (R 4.1.2)
scales 1.2.1 2022-08-20 [1] CRAN (R 4.1.2)
seqinr 4.2-30 2023-04-05 [1] CRAN (R 4.1.2)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.0)
shape 1.4.6 2021-05-19 [1] CRAN (R 4.1.0)
shiny 1.7.5 2023-08-12 [1] CRAN (R 4.1.2)
shinydashboard 0.7.2 2021-09-30 [1] CRAN (R 4.1.0)
stringi 1.7.12 2023-01-11 [1] CRAN (R 4.1.2)
stringr 1.5.0 2022-12-02 [1] CRAN (R 4.1.2)
SummarizedExperiment * 1.24.0 2021-10-26 [1] Bioconductor
survival 3.5-7 2023-08-14 [1] CRAN (R 4.1.2)
tibble 3.2.1 2023-03-20 [1] CRAN (R 4.1.2)
tidyr 1.3.0 2023-01-24 [1] CRAN (R 4.1.2)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.1.2)
tmvtnorm 1.5 2022-03-22 [1] CRAN (R 4.1.2)
tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.1.2)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.1.0)
usethis 2.2.2 2023-07-06 [1] CRAN (R 4.1.2)
utf8 1.2.3 2023-01-31 [1] CRAN (R 4.1.2)
vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.1.2)
vroom 1.6.3 2023-04-28 [1] CRAN (R 4.1.2)
vsn 3.62.0 2021-10-26 [1] Bioconductor
withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.2)
xfun 0.40 2023-08-09 [1] CRAN (R 4.1.2)
XICOR 0.4.1 2023-04-21 [1] CRAN (R 4.1.2)
XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.1.2)
xml2 1.3.5 2023-07-06 [1] CRAN (R 4.1.2)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
XVector 0.34.0 2021-10-26 [1] Bioconductor
yaml 2.3.7 2023-01-23 [1] CRAN (R 4.1.2)
zlibbioc 1.40.0 2021-10-26 [1] Bioconductor
zoo 1.8-12 2023-04-13 [1] CRAN (R 4.1.2)
[1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
R ── Package was removed from disk.
──────────────────────────────────────────────────────────────────────────────References
Yuan, Zuo-Fei, Simone Sidoli, Dylan M Marchione, Johayra Simithy, Kevin A Janssen, Mary R Szurgot, and Benjamin A Garcia. 2018. “EpiProfile 2.0: A Computational Platform for Processing Epi-Proteomics Mass Spectrometry Data.” Journal of Proteome Research 17 (7): 2533–41. https://doi.org/10.1021/acs.jproteome.8b00133.